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The magnetorotational instability (MRI) is thought to be a powerful source of 
turbulence and momentum transport in astrophysical accretion discs, but obtaining 
observational evidence of its operation is challenging. Recently, laboratory experi¬ 
ments of Taylor-Couette flow with externally imposed axial and azimuthal magnetic 
fields have revealed the kinematic and dynamic properties of the MRI close to the 
instability onset. While good agreement was found with linear stability analyses, lit¬ 
tle is known about the transition to turbulence and transport properties of the MRI. 

We here report on a numerical investigation of the MRI with an imposed azimuthal 
magnetic field. We show that the laminar Taylor-Couette flow becomes unstable 
to a wave rotating in the azimuthal direction and standing in the axial direction 
via a supercritical Hopf bifurcation. Subsequently, the flow features a catastrophic 
transition to spatio-temporal defects which is mediated by a subcritical subharmonic 
Hopf bifurcation. Our results are in qualitative agreement with the PROMISE ex¬ 
periment and dramatically extend their realizable parameter range. We find that as 
the Reynolds number increases defects accumulate and grow into turbulence, yet the 
momentum transport scales weakly. 
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I. INTRODUCTION 


The magnetorotational instability (MRI) is of great importance in astrophysics. First 
discovered by Velikhov [l| in 1959, it remained unnoticed until 1991 when Balbus and Hawley 
realised its application to accretion disc theory. Accretion discs are astrophysical systems 
that consist of ionised gas and dust orbiting a massive body. Planets and stars are formed 
from this initially dispersed matter. The physical mechanism of accretion is straightforward: 
a parcel of viscous fluid in the differentially rotating disc loses its angular momentum over 
time and falls onto the central object. To explain the astrophysically observed rates of 
accretion, however, one must assume a turbulent transport of angular momentum in the 
outward direction |3|. In so-called Keplerian discs the angular velocity prohle of gas follows 
the law 

Q ~ ( 1 ) 

which is hydrodynamically stable according to the Rayleigh criterion for rotating fluids 

d{r‘^VtY 


dr 


> 0 for stability. 


( 2 ) 


Ionized accretion discs, however, are necessarily magnetized and the MRI may still act in 
rotating flows, provided the angular velocity decreases with radius, which is true of Keplerian 
flows ([I]). 

The growth rates of the MRI and the parameter ranges in which it acts were determined 
in several linear analyses j^, [siJsl, but these do not provide information about the flow struc¬ 
ture and scaling of angular momentum transport after nonlinear saturation. In the last 
two decades there has been a great deal of numerical work concerned with the nonlinear 
properties of the MRI. Simulations are usually performed with the shearing sheet approxi¬ 
mation, which is a local model of an accretion disc with shear-periodic boundary conditions 
in the radial direction jb]. The main disadvantage of this model is the influence of boundary 
conditions on the geometry of the observed modes and transport scaling. In particular, the 
length of the computational box hxes the modes that appear and determines their nonlinear 
saturation. As it is not clear how the length should be selected, the interpretation of the 
transport scaling becomes quite involved |9|. 

These theoretical results and numerical simulations inspired physicists to realise the MRI 


in laboratory experiments. In 2001 Ji et al. 


10| and Rudiger and Zhang [m| independently 
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suggested the possibility of directly observing magnetorotational instabilities in a cylindrical 
vessel made of two co-axial and independently rotating cylinders containing a liquid metal 
alloy (see Fig. [T^). The standard form of the MRI (SMRI), which they proposed, emerges 
when a purely axial magnetic held is imposed, but this has not yet been achieved in ex¬ 
periments. The difficulty is that liquid metals have very small magnetic Prandtl numbers 
(e.g. Pm ~ 10“® for gallium alloys), leading in this case to very large Reynolds numbers 
{Re > 10^) necessary to observe the SMRI (see Table |T]for the dehnitions of Re and Pm). In 
fact, such high Reynolds numbers have never been achieved even for non-magnetic Taylor- 
Couette hows. A further difficulty of Taylor-Couette experiments in the quasi-Keplerian 
regime arises because of Ekman vortices that arise adjacent to the endplates. Unless a very 
specihc endplate arrangement is used, the Ekman vortices extend deep into the how and 
even at moderate Reynolds number the basic Couette how cannot be obtained experimen¬ 


tally 12|. The resulting velocity prohles are no longer quasi-Keplerian and 
instabilities render the how turbulent even in the absence of magnetic helds 


lydrodynamic 


13 


Hollerbach and Rudiger proposed instead a combination of axial and azimuthal mag¬ 
netic helds, giving rise to the helical MRI (HMRI) at much lower Re ~ 10^ for Hartmann 
numbers Ha 10 (see Table [I] for the dehnition of Ha). This was successfully observed 
15|, ll^ in the PROMISE facility (Potsdam-ROssendorf Magnetic Instability Experiment). 
Both the SMRI and HMRI consist of axisymmetric toroidal vortices, which are station¬ 
ary for the former but travel axially for the latter. Hollerbach et al. js] realised that al¬ 
though a purely azimuthal magnetic held does not yield any axisymmetric instabilities, 
non-axisymmetric modes can be destabilised. The resulting azimuthal MRI (AMRI) arises 
in Taylor-Couette how at Re 10^ for iPa ~ 10^ [st], and a recent upgrade of the PROMISE 
power supply made its experimental observation possible jl^ . In PROMISE the endplates 
are split into two parts and the inner (outer) one is attached to the inner (outer) cylinder. 
Although this conhguration is acceptable at Re < 3000, as studied experimentally, endplate 
ehects become dominant at larger Re. Because of this, and of the practical impossibility 


M- 


of generating a purely azimuthal magnetic held experimentally, it is chal 
identify the AMRI modes in the experimental data unambiguously (see 


enging to actually 


m- 


Despite this recent experimental progress in realising magnetorotational instabilities in 
the laboratory, little is known about their bifurcation scenario, transition to turbulence and 
transport properties as Re increases. In this work we address these points for the AMRI. We 
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perform direct numerical simulations of the coupled induction and Navier-Stokes equations 
using axially periodic boundary conditions, thereby avoiding undesired endplate effects and 
focusing on the features intrinsic to the AMRI. We hnd that the laminar quasi-Keplerian 
flow becomes unstable to a wave rotating in the azimuthal direction and standing in the 
axial direction. Subsequently, we identify a new bifurcation scenario giving rise to spatio- 
temporal defects via a subcritical subharmonic Hopf bifurcation. As the Reynolds number is 
further increased, the flow becomes turbulent and outward momentum transport is enhanced, 
albeit at a weak rate. The results are in good qualitative agreement with the PROMISE 
observations [l^ and substantially extend the parameter range explored experimentally. 


II. GOVERNING EQUATIONS 


We consider an incompressible viscous liquid metal that is sheared between two inde¬ 
pendently rotating cylinders of radii (inner) and Tq (outer). The angular velocity of the 
cylinders are Qj and Qq, respectively, and an external azimuthal magnetic held {ri/r)Bo, 
where r is the radial coordinate, is imposed. The relevant huid properties are the electrical 
conductivity a, the kinematic viscosity u , the density p and the magnetic diffusivity 77 . The 
velocity held u is determined by the Navier-Stokes equations (|3]), whereas the magnetic held 
B is determined by the induction equation (jlj), which represents a combination of the laws 
of Ampere, Faraday, and Ohm. The equations were rendered dimensionless by using the 
gap between cylinders d = Tq — Vi for length, d^/z/ for time, and Bq for the magnetic held. 
In dimensionless form they read 


(dt + v ■ V)v = —Vp V^v -h —(V X B) X B, 

Bm 

(ai-^V2)B = Vx(vxB), 

Bm 


( 3 ) 


( 4 ) 


together with V ■ v = V ■ B = 0. Here p is the pressure, Ha the Hartmann number, and Bm 
the magnetic Prandtl number. The dimensionless parameters of the system are specihed in 
Table m Following the PROMISE experiment {^, we use a magnetic Prandtl number of 
Bm = 1.4 • 10“® (corresponding to the alloy Ga^'^, a radius-ratio of <5 = 0.5 and 
a rotation-ratio of p = 0.26. This places the velocity prohle in the quasi-Keplerian regime, 
for which the angular velocity decreases radially, whereas the angular momentum increases, 
i.e. (5^ < /i < 1. 
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(a) (b) 

B 




FIG. 1. (a) Schematic of the Taylor-Couette geometry with azimuthal magnetic field. 

A liquid metal is confined between two coaxial cylinders of radii and Tq, which can rotate 
independently at angular velocities and Hq. In this work the rotation rate is fixed at /r = 
Qo/^i = 0.26 and an azimuthal magnetic field of the form (rj/r)i?o is imposed, where r is the 
radial coordinate. ( 6 ) Instability region of the AMRI. Blue circles correspond to the points at 
which our simulations were conducted, the red dashed line to the curve of maximum growth rate 
and the green solid line is a fit to the latter of the form Ha = aRe^, where a = 0.71 and b = 1.55. 

A. Boundary conditions 


We employ cylindrical coordinates (r, 0, z) G [r,, To] x [0, 27r] x [0, L^], for which the no-slip 
velocity boundary conditions at the cylinders read 


0, z) = Re, u^ir^, 4>, z) = fiRe. 


( 5 ) 


Periodicity in the axial direction is imposed with basic length L^. The background circular 
Couette flow V = l/(r)e 0 is a solution to the equations and boundary conditions given by 


V{r) = 


1 


1 + 5 




-(i?e — fi Re)- 


The magnetic field is also assumed to be periodic in the axial direction, 
direction the boundary condition depends on the material of the cylinders. 


( 6 ) 

In the radial 
Typically two 
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TABLE I. Dimensionless parameters of the magnetohydrodynamic Taylor-Couette 
problem. 


Abbrev. 

Parameter 

Definition 

Range 

5 

Radius ratio 

n/vo 

0.5 

a 

Axial wavenumber (geometrical parameter) 

2-nlL^ 

0.5—4.5 

h 

Angular velocity ratio 


0.26 

Pm 

Magnetic Prandtl number 

vjri 

1.4 • 10"® 

Re 

Reynolds numbers of inner cylinder 

^iTidjv 

1480—9333 

Ha 

Hartmann number 


90—457 


idealized cases are considered in the MRI problem; insulating and conducting cylinders. 
These lead to slightly different results, as theoretically demonstrated by Chandrasekhar 18 |. 
However, the difference is not great, and here we will consider only the case of insulating 
boundaries. Assuming the Fourier expansion for each variable 


A = ^ ^ Afc,m(r)exp[i(Q;fc^ + m0)], (7) 

\k\<K \m.\<M 

one obtains the following boundary conditions for B: 

For case k = m = 0: 

Be = Bz = 0 . ( 8 ) 


case k = 0, m ^ 0: 


Br ± iBg = 0, = 0 (+ on r*, - on Tq). 


( 9 ) 


case k ^ 0: 


Om{c(kn) r 


( 10 ) 


where Bm{x) denotes the modified Bessel 
and = d^Bm- See Willis and Barenghi 


'unction Im{x) for r = ri and Km{x) for r 


19| for a detailed derivation. 


= r 


O? 
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B. Brief remarks on symmetries of rotating magnetohydrodynamic flows 


The basic circular Couette flow ([6]) has SO(2) x 0(2) symmetry, where SO(2) represents 
the rotational symmetry in the azimuthal direction. In the axial direction the group 0(2) 
may be written as 0(2)= Z 2 xSO( 2 ), where Z 2 is a reflection (up-down symmetry) and 
SO(2) the translational symmetry in the 2 ; direction. The presence of purely axial or purely 
azimuthal imposed magnetic held does not change the symmetry group of the system. Hence 
if the primary instability is a Hopf bifurcation the resulting states can be either standing or 
traveling waves in the axial direction j^|. By contrast, a combined helical magnetic held 
breaks the reflection symmetry and only traveling waves (TW) can be observed 2l|]. Finally, 
if the bifurcating solution is non-axisymmetric, as in the AMRI, this will generically be a 
rotating wave in the azimuthal direction. 


III. NUMERICAL METHOD 


In the numerical simulations only the deviation from the basic flow u = v—V is computed. 
Its governing equations read 


(a* - V^) u = N - Vp, V ■ u = 0, 


( 11 ) 


which are supplemented with homogeneous boundary conditions u = 0. Here N stands for 
the nonlinear term in the Navier-Stokes equations ([3]), which contains the advective terms 
and the Lorentz force: 

Q 


N = u X (V X u) - (V • V)u - (u • V)V 


Pm 


(V X B) X B 


( 12 ) 


Q 


= u X (V X u) — (l//r)c?^u -h {2Vlr)u^er — Ur{l + dr)Ve^ + —— (V x B) x B. 


Spatial discretisation is accomplished via the Fourier expansion in the azimuthal and axial 
directions o. and as the variables are real, their Fourier coefficients satisfy the property 
Ak^m = A*_^ where A* denotes the complex conjugate. 

The pseudospectral Fourier method is the most efficient choice for periodic boundary 
conditions. Because of their great accuracy at low resolutions, spectral methods have also 
been used to discretise the hydromagnetic Taylor-Couette flow problem in the radial di¬ 


rection 


19 


22| . In these works the Chebyshev collocation method was chosen because of 









its simplicity. Nevertheless its computational and storage costs scale as where N 

is the number of radial points, making computations at large Reynolds number imprac¬ 
tical. Petrov-Galerkin formulations reduce the cost to 0{N) and have been also used in 


hydrodynamic Taylor-Couette flows 


m 


24l | . However, the treatment of the radial boundary 


conditions becomes very cumbersome for the magnetic held. We here use the hnite-difference 
method in the radial direction. Radial derivatives are calculated using 9-point stencils to 
(9 — order, where j is the order of the derivative. This results in banded matrices 
with associated 0{N) cost, while providing excellent accuracy. Ref. |25| provides a more 
thorough discussion of these computational issues and the accuracy of the hnite-difference 
method applied to hydrodynamic Taylor-Couette how. 

A second-order scheme is applied to the time discretisation = qAt, based on the 
implicit Crank-Nicolson method. Applying this method to the Navier-Stokes equations, we 
hnd 

(1/At-= (1/At+ V2 )u‘' + N'?+^-Vp, vV = V-N''+T (13) 

where is an estimate for the nonlinear terms (Euler predictor, Crank-Nicolson correc¬ 

tor). In cylindrical coordinates the r- and 0-components of the Laplacian operator couples 
the two components, and for a Fourier decomposition they are complex operators. Program¬ 
ming complexity and computational cost of inversion for and can be reduced, 
however, by considering 


U± = Ur A iu^, 


i.e. Ur = lr{u+ + u-), - U-), 


(14) 


where the ± are taken respectively. The equations governing these components separate 
and the Laplacian operator is now real im). 


(0i-V^)n± = A±-(Vp)±, 


= V" 


4 ± Xd.. 


(15) 


A. Influence-matrix method 


The natural approach to solving (lT3il is to invert for p and then for All the bound¬ 

ary conditions are on however, there are none on p, and it is well known that primitive 


variable formulations are su 
ditions are enforced on p 


□ject to loss of temporal order if inappropriate boundary con- 


26| . For the magnetic held, there appear at hrst sight to be too 
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few boundary conditions, and further, the components of B are coupled in the boundary 
condition. It is shown here that the influence-matrix method resolves these issues, the ap¬ 
propriate boundary conditions can be satished to machine precision, and temporal order is 
retained. We show hrst how the method is applied for time integration of the velocity held, 
similar to j^. An analogous approach is applied to the magnetic held. 


1. Method for the velocity field 


We write the time-discretised Navier-Stokes equations (IT^ in the form 


= Fu-? + N‘'+5 - Vp, 
= V-(Fu^ + N'^+^), 


(16) 


where q denotes time tq. This form is sixth order in r for and second order for p, without 
the solenoidal condition explicitly imposed. In principle this system should be inverted 
simultaneously for p and with boundary conditions = 0 and V ■ = 0. In 

practice it would be preferable to invert for p hrst then for but the boundary conditions 
do not involve p directly. Note that the Yvfl term has been included in the right-hand side 
of the pressure-Poisson equation, so that it corresponds precisely to the divergence of the 
equation for This ensures that any non-zero divergence in the initial condition is 

projected out after a single time-step. We split the system flT6|) into the ‘bulk’ solution, 

Xu = yu"? -|- — Vp, 

v^p = v-(yu'^ + N‘?+^), 

with boundary conditions u = 0 and d^p = 0, and the auxiliary systems 


(17) 


A'u' = -Vp', 

VV = 0. 


(18) 


with two sets of boundary conditions, u' = 0 with drp' = {0,1} and {1,0} on r = {rj,ro}, 
and 

{ Xu' = 0, (19) 

with boundary conditions = {0,1}, u'_ = {0,1}, u'^ = {0, i} and similarly their reversed 
versions on {r*. To}- These dashed functions may be precomputed, and will be used to correct 
the approximate boundary conditions used to calculate u and p = 0. 
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The system ffTSD provides, with the two boundary conditions, two linearly independent 
functions u'- that may be added to u without altering the right-hand side in flTTl) . Similarly 
the system flT^ provides a further six functions. The superposition 

8 

= u + ^ Oj u' (20) 

i=i 

may be formed in order to satisfy the eight original boundary conditions, = 0 and 

V ■ = 0 on r = {rj,ro}. Substituting (120|1 into the boundary conditions, they may be 

written 

2la=-g(u), (21) 

where A = 74(g(u')) is an 8x8 matrix. The appropriate coefficients required to satisfy the 
boundary conditions are thereby recovered from the solution of this small system for a. The 
error in the boundary conditions using the influence-matrix technique is at the level 

of the machine precision. 

The auxiliary functions u'j{r), the matrix A and its inverse may all be precomputed, and 
the boundary conditions for u' have been chosen so that m'_i_ are purely real, u'^ is purely 
imaginary, and A is real. For each timestep, this application of the influence matrix technique 
requires only evaluation of the deviation from the boundary condition, multiplication by an 
8x8 real matrix, and the addition of four functions to each component of u, each either purely 
real or purely imaginary. Compared to the evaluation of nonlinear terms, the computational 
overhead is negligible. 


2. Method for the magnetic field 

Consider the induction equation (|4]) time-stepped without V ■ B = 0 enforced. Evolution 
of "0 = V ■ B is then governed by the divergence of (|4|) . 


dtif = 

Fm 


( 22 ) 


In addition to the boundary conditions derived in 


19| . the condition '0 = 0 must be sat- 
ished on the boundary to avoid introduction of divergence into the domain. Then (jl]) has 
the appropriate number of boundary conditions, and 0 = V ■ B should remain zero for a 
solenoidal initial condition. 
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To prevent accumulation of divergence from artificial internal sources, i.e. discretisation 


error, it is commonplace to introduce an artificial pressure If 


The discretised system 


is then as in flTHl) where one reads B for u and If for p. The boundary condition for fl is 
any choice such that, when one computes If for a given B'^, it is found to be constant when 
V • B”? = 0 The choice = 0 is applied here. When comparing with the problem for 
the velocity, here the difficulty is not the coupling of the boundary condition for B with If, 
but between the components of B at an insulating boundary. 

Here the system is split as in ffT7|) for the ‘bulk’ solution, with approximate boundary 
condition B = B*^ on {rj,ro}. This is then corrected precisely via the influence matrix 
requiring only the simple auxiliary system 


XB' = 0, 


(23) 


with boundary conditions B'^ = {0,1} or {1, 0} and = {0, i} or {i, 0} on {rj. To}- 

Problem fl23|) separates for the three components, which, with the two boundary condition 
options for each, provides six functions Bj{r). The correction is then 

6 

B'^+i = B + ^a,B'. (24) 

Let gj{B) = 0 denote the insulating boundary conditions and solenoidal condition evaluated 
at Tj and r^. Substituting fl2T)) into the boundary conditions, they may be written Ha = 
-g(B), where A = H(g(B')) is a 6x6 matrix. The appropriate coefficients required to 
satisfy the boundary conditions are recovered from solution of this small system for a. 

Again, the auxiliary functions B'(r), the matrix A and its inverse may be precomputed, 
and the boundary conditions for B' have been chosen so that B'^ are purely real, i?' is purely 
imaginary, and A is real. At the end of the timestep, the solution is solenoidal and satisfies 
the boundary conditions to machine precision. 


B. Implementation notes and parallelization 


The Taylor-Couette flow code was written in FortranQO. Nonlinear terms are evalu¬ 


ated using the pseudo-spectral method and are de-a! 


transforms are performed wit 


1 the FFTW3 library 


iased using the 3/2 rule. The Fourier 


30| and matrix and vector operations 


are performed with BLAS 3l|. Each predictor-corrector iteration involves the solution of 
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banded linear systems with forward-backward substitution using banded LU-factorizations 
that are precomputed prior to time-stepping. These operations are performed with LAPACK 
32| . The code was parallelized so that data is split over the Fourier harmonics for the linear 
parts of the code: evaluating curls, gradients and matrix inversions for the time-stepping 
(these linear operations do not couple modes). Here all radial points for a particular mode 
are located on the same processor; separate modes may be located on separate processors. 
Data is split radially when calculating Fourier transforms and when evaluating products in 
real space (nonlinear term of the equations). The bulk of communication between processors 
occurs during the data transposes. 


C. Numerical validation 

The code was validated against several published linear stability results, as well as three- 
dimensional nonlinear simulations of the coupled induction and Navier-Stokes equations. 
We tested the inductionless limit Pm = 0 and hnite Pm, obtaining excellent agreement in 
all cases. 


i. Linear stability of Couette flow subject to magnetic fields 


Linear instabilities were detected in the calculations by monitoring the kinetic energy of 
the deviation from circular Couette flow after introduction of a small disturbance. In the 
linear regime we write 

u' ~ exp(At -|- \[kz + m0]), E ~ \u'\^ ~ exp(2(Tt), 


where A = a -I- icu is a complex number; the imaginary part uj is the oscillation frequency and 
the real part a the growth rate of the dominant perturbation. The latter is readily extracted 
from the relationship log(P) ~ 2at. 

We hrst reproduced the classical results of Roberts (33| , who considered the inductionless 
limit Pm = 0 for narrow gap geometry 5 = 0.95 and stationary outer cylinder. For a 
Hartmann number of Pa = 5.477 he obtained a critical Reynolds number of Rcc = 281.05 
with associated critical axial wavenumber of Oc = 2.69 and m = 0. In our simulations we 
hxed a = ac and obtained Rcc = 281.055 using N = 33 radial points. For hnite magnetic 
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Prandtl number Pm = 1 we reproduced the results of Willis and Barenghi 3^ for wide gap 
5 = 0.5 and stationary outer cylinder. For Ha = 39, and a = 2.4 and m = 0 they found 
Rcc = 60.5, which is in good agreement with our result {Rcc = 60.3). In order to test the 
azimuthal magnetic field we reproduced recent results of Hollerbach et al. js] for the AMRI. 
For example at Pm = 0, Ha = 316, Re = 1000, 6 = 0.5 and /i = 0.26, they obtained 
a = —78.6 for wavenumbers a = 7.17 and m = 1, which is in very good agreement with our 
value a = —78.7. 


2. Nonlinear simulations 


Willis and Barenghi 35|] explored dynamo action in Taylor-Couette flow. They first 
solved the Navier-Stokes equations in the absence of magnetic held and subsequently applied 
a small magnetic disturbance to test whether it grew into a dynamo. In the axisymmetric 
Taylor-vortex regime axisymmetric magnetic helds were found to decay, in accordance with 
Cowling’s anti-dynamo theorem. Non-axisymmetric magnetic helds may be excited, how¬ 
ever, and for Re = 136.4, 6 = 0.5, Pm = 2, a = 1.57 and stationary outer cylinder they 
observed that the magnetic disturbance grows for m = 1 {aB,m=i ~ 0 . 2 , leading to dynamo 
action), whereas it decays for m = 2 {aB,m =2 ~ —1.4). We reproduced this setting using 
N = 41, K = 16 and M = 12 and obtained (TB,m=i ~ 0.16 and aB,m,=i ~ —1.42, in good 


agreement with 


35| 


Finally, we compared results of the axisymmetric HMRI (helical held B = Ro(Gz + 7e<^)) 


obtained with the spectral code of Hollerbach 
torque at the cylinders 


22 | . A typical diagnostic quantity is the 


o d 

G ~ —27rr^ — 
or 


r 


2'Kr^ 


o 

— - OrU^ 

r 


( 26 ) 


The laminar how torque will be used as a scale, so that the dimensionless ratio G/Giaminar 
measures the intensity of angular momentum transfer relative to laminar how. We choose the 
parameters Re = 300, Ha = 10, <5 = 0.5, 7 = 2, a = 0.314, which are well into the nonlinear 
regime. After nonlinear saturation the dimensionless torque on the cylinders obtained with 
our code for A^ = 81 and K = 192 was G/Giaminar = 1.4122, which is in excellent agreement 
with the code of Hollerbach (G/Giaminar = 1.4123). 
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IV. PRIMARY INSTABILITY: STANDING WAVES 


In the experiments of Seilmayer et al. 17| the AMRI was explored near the onset of 
instability for two different Reynolds nnmbers Re = 1480 and 2960, and Hartmann numbers 
in the range Ha G [0,160]. The experiments have an aspect-ratio of 10. Here we selected a 
periodic domain of length = 12.6, and initialised the simulations by disturbing all Fourier 
modes with the same amplitude, thus allowing the axial wavenumber to be naturally selected. 
Because of the symmetries (see §II Bll . two different Hopf-bifurcation scenarios are possible 


2l| . In the hrst one, the z-reflection symmetry is broken and depending on the initial 


conditions either upward traveling waves (with /c > 0 modes) or downward traveling waves 
(with /c < 0 modes) may be observed. In the second scenario, the ^-reflection symmetry is 
preserved and a standing wave emerges. This is a combination of upward and downward 
traveling waves for which positive and negative k modes are in phase and have exactly the 
same amplitude. In both scenarios waves rotate in the azimuthal direction. 


We found that at Re = 1480 the circular Couette flow ([6]) becomes unstable at Hac = 107. 
The emerging pattern is a standing wave (SW) with dominant mode {k,m) = (±8,1), so 
that 8 pairs of vortices £t in the domain. Figure |2 ]d shows the square of the amplitude of 
the complex Fourier coefficient Ag i for increasing Ha. As expected in a Hopf bifurcation. 
Aim oc Ha — Hac near the onset of instability, and this relationship holds up to FTa ~ 
112. The vortex arrangement of the standing wave at Ha = 150 is shown in the flow 
snapshot of hgure |2^. In this case the mode k = 9 was naturally selected. Thus the 
dominant axial wavenumber depends on Ha because of the Eckhaus instability, as also 


observed in hydrodynamic Taylor-Couette flow 


36| . The torque changes respectively with 


axial wavenumber (black curve on the Fig. [3^), so at the same parameter value states 
with different wavenumber and torque can be realised depending on the initial conditions. 
Further increasing Ha the instability is gradually damped until it disappears at Ha ~ 175. 
Over the whole Hartmann range the additional torque due to the SW never exceeds 1% of 
the laminar flow (see hgure [3^), indicating very weak transport of angular momentum. The 
maximum in torque correlates well with the maximum growth rate from the linear stability 
analysis shown in Fig. [3 ]d. 
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FIG. 2. Primary instability: a standing wave arises through a supercritical Hopf bifur¬ 
cation. (a) Standing wave with k = 9 (corresponding to 9 vortex-pairs in the axial direction) 
at Re = 1480, Ha = 150 and = 12.6 (long domain). From left to right: isosurfaces of axial 
velocity Vz = ±0.005 (normalized with the velocity of the inner cylinder rijr*), contours of axial 
and radial velocity. The aspect ratio of the colormaps has been stretched by a factor of 0.6. (b) 
Onset of instability at Re = 1480. The critical Hartmann number is Huc ~ 107 with critical 
axial mode k = 8. The square of the amplitude of the Fourier coefficient ^ depends linearly on 
Ha — Hac close to the critical point as expected in a Hopf bifurcation. The coefficient H_8^i has 
the same amplitude as confirming that the axial reflection symmetry is preserved (standing 
wave). 

V. ONSET OF SPATIO-TEMPORAL CHAOS 

At Re = 2960 a Hopf bifurcation occurs at Hqc = 120 and the emerging SW remains 
stable until Ha = 160. Increasing Ha beyond this point a catastrophic transition to spatio- 
temporal chaos is observed: the vortex structure is damaged and the up-down symmetry 
is broken (Fig. 0]). Between Ha = 130 and 160 there is a hysteresis region in which both 
SW and spatio-temporal chaos (defects) are locally stable (see Fig. [3^). In this i7a-range, 
if the initial condition is a SW from another run with slightly different Ha, this remains 
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(a) (b) 



FIG. 3. Onset of spatio-temporal chaos, (a) Dimensionless torque for AMRI versus Ha 
for Re = 1480 and Re = 2960. Eckhaus instability at Re = 1480: the branches of the black 
curve belong to different axial wavenumbers (k=8, 9 and 10) of the standing wave. Bistability 
at Re = 2960: in the yellow-shaded region standing waves (green) and defects (red) coexist; 
between them there is an unstable branch or edge state (blue), (b) Perturbation growth rates 
a (normalised with llj) as a function of Ha for Re = 1480 and Re = 2960. Positive values of a 
correspond to instability. 

stable. However, when starting for example from a randomly disturbed Couette profile the 
flow evolves directly to defects. 

This catastrophic transition suggests a subcritical bifurcation. We investigated this hy¬ 
pothesis by computing the unstable branch separating defects and SW. For this purpose we 
combined time-stepping with a bisection strategy as follows. If the SW is slightly disturbed, 
then the flow should rapidly converge to the SW because it is locally stable. The same ap¬ 
plies to defects. For intermediate initial conditions the flow should take a long time before 
asymptotically reaching either the SW or the defects. Such initial conditions were generated 
here by performing a linear combination between two selected flow snapshots of SW and 
defects. This combination was parametrised with a variable /S, for which /3 = 0 corresponds 
to SW and /3 = 1 to defects. With the bisection procedure, rehning (5 results in an initial 
condition successively closer to the manifold (or edge) delimiting the two basin boundaries. 
The edge is comprised of those initial conditions that tend neither to defects nor to SW, 
and the attractor in this manifold is referred to as an edge state . 
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(a) (b) 



FIG. 4. (a) Defects at Re = 2960,77o = 190 and = 12.6 (long domain). From left to right: 
isosurfaces of axial velocity {vz = ±0.01 [riirj]), contours of axial and radial velocity, {h) Onset 
of turbulence. Isosurfaces of axial velocity {vz = ±0.0125 [fljr*]), contours of axial and radial 
velocity. Re = 4000, Ha = 264 and Lz = 12.6. The aspect ratio of the colormaps has been 
stretched by a factor of 0.6. 

Figure [5] shows that as initial conditions are taken closer to the edge, the temporal 
dynamics become simple as the edge state is approached. At Ha = 155, which is very close 
to the destabilisation of the SW, the dynamics appear to exhibit a damped oscillation (see 
Fig. [5b). Unfortunately, it is difficult to establish whether the oscillation hnally decays or 
saturates at a tiny amplitude, as expected close to the bifurcation point. At Ha = 140, 
however, which is further from the bifurcation point, the oscillation saturates at non-zero 
amplitude (see Fig. |5b). This suggests that the edge state is a relative periodic orbit (or 
modulated wave) emerging at a subcritical Hopf bifurcation of the SW. Despite this simple 
temporal behaviour, the spatial structure of the edge state is complicated (see Fig. |6b-c). It 
consists of a long-wave (subharmonic) modulation of the axially periodic pattern of the SW, 
which can be seen as a precursor to defects (compare Fig. [6b-b to Fig. IHi). We expect that 
as Ha is further reduced the edge state suffers a bifurcation cascade and becomes chaotic. 
This should continuously connect to defects and stabilise at a turning point for Ha > 130, 
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(a) (b) 




FIG. 5. Edge tracking procedure at Re = 2960. (a) Ha = 140 and (b) Ha = 155. Green 
curves evolve toward the standing wave state and red curves toward defects, although all of them 
start very close to the edge state. Oscillations at Ha = 155, which is close to the destabilisation 
point of the standing wave, appear to decay, while at Ha = 140 they saturate. Time is normalised 
using the inner cylinder rotation frequency 1/lli, i.e. t = t ■ Re. 

which is the lowest Ha for which defects remain stable. 

VI. TURBULENT TRANSPORT OF MOMENTUM 

As the Reynolds number is further increased, defects are expected to grow gradually into 
turbulence. Although it would be very interesting to perform a two-parameter study of the 
dynamics m Ha and Re, this is computationally expensive and beyond the scope of the 
current work. We here chose to follow a parameter path of the form 

Ha = aRe^, (26) 

with a = 0.71 and b = 1.55. This path is shown as a solid green line in Fig. [T}3 and provides 
a very good approximation to the curve of maximum growth rate of the linear stability 
analysis (red dashed line). It goes deep into the instability region and so we expect the 
instability to fully develop as Re increases with Ha subject to (126|1 . 

At Re = 4000 the vortices are small at the inner cylinder and remain quite large at the 
outer cylinder (Fig. |4 )d), and at Re = 9333 this tendency develops into rapidly drifting small 
vortices at the inner cylinder and slow large vortices at the outer cylinder (Fig. [7^). There 
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FIG. 6. Edge state at Re = 2960 and = 12.6. (a) Close to the bifurcation point, 

Ha = 155: isosurfaces of axial velocity = ±0.0135 [fljrj]. (b) Far from the bifurcation 
point, Ha = 140: isosurfaces Vz = ±0.012 [fliri]. (c) Contours of axial and radial velocity, 
Ha = 140. The aspect ratio of the colormaps has been stretched by a factor of 0.6. The edge state 
consists of a long-wave modulation of the standing wave. 


is no preferred direction in the system; vortices can travel up or down, both at the inner 
and outer cylinders. 

The qualitative difference between standing wave, defects and turbulent flow is apparent 
in time series of the radial velocity Vr taken at the mid-gap between the cylinders (r, 0, z) = 
(1.5, 0, 0). FigurelHk shows that the radial velocity of the standing wave oscillates periodically 
around zero. The edge state features a slow temporal frequency modulating the oscillation of 
the SW (Fig. [8 )d). For defects at Re = 2960 the time series is mildly chaotic. As Re increases 
toward turbulence the velocity pulsates in a very chaotic manner (Fig. (HJi). However, the 
main frequency associated with the AMRI can still be discerned. By comparing all panels 
it becomes apparent that this frequency scales with the rotation-rate of the inner cylinder. 
This is consistent with the linear stability analysis of 8|], and with the studies 38|, l39|, where 


it is shown that in the low Pm limit the AMRI is an inertial wave. 


The transfer-rate of angular momentum between the cylinders is important for accretion- 
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(a) (b) 




FIG. 7. (a) Turbulent flow in a short domain at Re = 9333, Ha = 456.7 and = 1.4. Axial 
velocity isosurfaces Vz = ±0.011 [fliri], contours of axial and radial velocity, {h) Dimensionless 
torque as a function of Re along the parameter path Ha = aRe^, with a = 0.71 and b = 1.55. 
This path is shown as a green line in Fig. [TJ). 


disc modelling. We checked the torque scaling for increasing Re numbers (see Fig. [7]o) 
along the parameter path (|26|) . The dimensionless torque increases with Re according to 
the scaling law 

G ~ (27) 


which is surprisingly low compared to hydrodynamic experiments in the Rayleigh unstable 
regime [^. We believe that this torque scaling comes close to being an upper bound for 
the torque scaling of the AMRI in the Re range studied here because the maximum in 
torque correlates well with maximum growth rates at low Re (see Fig. ED. However, we 
must caution that in hydrodynamic Taylor-Couette flow different maxima of the torque 


have been observed as a function of the relative rotation of the cylinders 


^|. At large Re 


these are not correlated to the maximum growth rate of the primary instability. Similar 
phenomena may occur for the AMRI. 


VII. DISCUSSION 


We showed that the AMRI in Taylor-Couette flow manifests itself as a wave rotating in 
the azimuthal direction and standing in the axial direction, thereby preser ving the reflection 


symmetry in the latter. In order to compare to experimental observations 


17j we computed 
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FIG. 8. Transition to turbnlence. Evolution of radial velocity perturbation Ur at the point 
{r,(p,z) = (1.5,0,0) with time (a) at Re = 1480 and Ha = 150 (standing wave); (6) - the same 
at Re = 2960 and Ha = 140 (edge state); (c) - at Re = 2960 and Ha = 190 (defects); (d) - at 
Re = 9333 and Ha = 456.7 (turbulence). Time is scaled using the inner cylinder rotation frequency 
l/flj, i.e. t = t ■ Re; velocities are normalised with the velocity of the inner cylinder 11 jr*. 


the angular drift frequency of the wave. This is shown in £gure|9] after being normalised with 
the rotation frequency of the inner cylinder. The wave rotates at approximately the outer 
cylinder frequency (dashed line) and slows down as the Hartmann number increases, which is 
in qualitative agreement with the experimental data. Note, however, that in the experiment 
two frequencies are simultaneously measured, corresponding to the up- and down-traveling 
spiral waves, respectively. Although in the standing wave the two frequencies are identical, 
in the experiment the asymmetric wiring creates and components of magnetic held 
which break the rehectional symmetry. As a result, up and down spirals travel with different 
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FIG. 9. Comparison to PROMISE experiment: angular drift frequencies of the waves at 
(a) Re 1480 and {h) Re 2960. Blue and red lines correspond to experimental results, black to our 
nonlinear simulations; the green line denotes outer cylinder rotation = 0.26. The waves 

rotate at approximately the outer cylinder frequency and slow down with increasing Ha. 


frequencies, similar to co-rotatin g T aylor-Couette flow in which the reflection symmetry is 


broken by an imposed axial flow 


42| . Another difference is that in the experiment the flow 


becomes unstable at lower Ha, which may be explained by the different boundary conditions 
in the experiment from our simulations. In the experiment copper cylinders are used, so 
perfectly conducting walls would be a closer boundary condition for the magnetic held. 
More signihcantly, in the experiments the cylinders are of hnite length, so to reproduce their 
results exactly a no-slip condition on end-plates should be used. We have applied periodic 
boundary conditions in the axial direction, which more accurately model the accretion disc 
problem and allow us to compute high Reynolds number hows more efficiently. 


As Re increases, a catastrophic transition to spatio-temporal chaos occurs directly from 
the SW. In a range of parameters SW and chaos are both locally stable and can be realised 
depending on the initial conditions. We have shown that the hrst step in this transition 
process is a subcritical Hopf bifurcation giving rise to an unstable relative periodic orbit, 
which has been computed using an analogue of the edge-tracking algorithm introduced by 
Skufca et al. in shear hows. This unstable relative periodic orbit consists of a long¬ 
wave modulation of the axially periodic pattern of the standing wave and destroys the 
homogeneity of the vortical pattern. It can thus be seen as a temporally simple defect 
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precursor of the ensuing spatio-temporal chaos. Because of the computational cost we could 
not track further instabilities on the unstable branch, which we speculate result in chaotic 
flow before the dynamics stabilize at a turning point {Ha = 130 at Re = 2960). After the 
turning point defects are stable and can be computed simply by time-stepping. 

We believe that such long-wave instabilities are ubiquitous in fluid flows. In linearly 
stable shear flows, such instabilities of traveling waves were found to be responsible for 


spatial localisation 
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44| . In fact, in pipe flow the ensuing localised solutions, which 


are also relative periodic orbits, suffer a bifurcation cascade leading to chaos 45(]. One 
difference is that in pipe flow the traveling waves are disconnected from laminar flow, where 
the standing wave of the AMRI is connected to the circular Couette flow. 

Our simulations were performed with a powerful spectral DNS method, which we have 
developed and validated with published results to excellent agreement. The method allowed 
us to compute flows up to Re = 10"^. As Re increases, defects accumulate and the flow 
becomes gradually turbulent. Although we found that the AMRI exhibits a weak scaling 
of angular momentum transport, with G oc Re^'^^, we expect that larger magnetic Prandtl 
numbers Pm (realistic for accretion discs) may result in a stronger scaling. Astrophysically 
important issues such as the precise angular momentum transport scalings obtained for 
different values of Pm, and also for different choices of imposed held (e.g. SMRI, HMRI, 
AMRI) will be the subject of future investigations. 
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